﻿using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Text;
using System.Threading;
using System.Xml.Serialization;

using MathNet.Numerics.Interpolation;
using MathNet.Numerics;
using MathNet.Numerics.Transformations;


namespace IdeaNMR.Shared
{
    public class BenchmarksBuilder
    {
        public volatile string ExecutionStatus;

        public volatile bool Terminated;

        private Thread Thread { get; set; }

        // ALISTAIR: COMMENT ADDED 20/05/2011
        // Declare array of BenchmarkSetup objects
        // where each BenchmarkSetup represents an 
        // indirect dimension of the experiment
        public BenchmarkSetup[] InputData { get; set; }
        
        // a list of analytical results, one set of results for each indirect dimension
        public List<BenchmarkResult> Results { get; private set; }

        // parameter specifying the rotational correlation time (tauC)
        // for which the user desires analytical data
        public double userTauC = 0;

        // ALISTAIR 06/03/2012: list of simulated tauC values, used for interpolation across tauC
        public List<double> tauCsForInterpolation;

        // ALISTAIR 06/03/2012: two dimensional list of points used to interpolate across TauC values
        List<List<Point>> timeDomainInterpolatedData;

        // Main method for processing each set of simulated FID data, applying the Fourier transform
        // and analysing the resulting spectra.
        // This method is called for each value of tauC for which Spinach data is supplied in XML format,
        // for each indirect dimension.
        // It is also called for the interpolated data at the user specified tauC, for each indirect dimension
        private BenchmarkResult CalculateBenchmarks(BenchmarkSetup bS, SimulatedData simData)
        {
            // ALISTAIR: COMMENT ADDED 24/01/2011
            // Read in transfer function as an array of points from BenchmarkSetup XML tag
            // (via BenchmarkSetup class) and order by time vaule ('T')
            Point[] TransferFunction = simData.GetTransferFunction().OrderBy(tf => tf.T).ToArray();
            
            // Diagnostic output for debugging
            /*
            String diagnosticOutput = "Points from Transfer Function:\n";
            diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
            diagnosticOutput += "Datapoint,Time,Real Value, Imag Value\n";
            for (int k = 0; k < TransferFunction.Length; k++)
            {
               diagnosticOutput += (k+1) +  "," + TransferFunction.ElementAt(k).T.ToString()
                                + "," + TransferFunction.ElementAt(k).V.Real.ToString()
                                 +"," + TransferFunction.ElementAt(k).V.Imaginary.ToString() + "\n";
            }
            System.Diagnostics.Debug.WriteLine(diagnosticOutput);
            */

            // Create an interpolation function
            // ALISTAIR: COMMENT ADDED 23/05/2011
            // Using natural cubic spline interpolation
            // The first parameter is the list of times in the FID
            // The second parameter is the list of signals in the FID
            // Note: Since the interpolation method can not handle complex numbers,
            // the real and imaginary components of each complex number are interpolated independently
            IInterpolationMethod realDataInterpolation = Interpolation.CreateNaturalCubicSpline(
                                TransferFunction.Select(tt => tt.T).ToList(),
                                TransferFunction.Select(tt => (double)tt.V.Real).ToList());
            IInterpolationMethod imagDataInterpolation = Interpolation.CreateNaturalCubicSpline(
                                TransferFunction.Select(tt => tt.T).ToList(),
                                TransferFunction.Select(tt => (double)tt.V.Imaginary).ToList());

            // ALISTAIR: COMMENT ADDED 20/05/2011
            // Interpolation of the transfer function resulting in
            // the number of points specified by the TimeDomain parameter.
            // ALISTAIR: MODIFICATION 09/03/2012
            // Change implementation of interpolated points to use calculated increment time (1/spectral width)
            // rather than being based on Tmax (but need to restrict Tmax as an upper limit 
            // for the increments that can be interpolated from the data available from Spinach)

            // list of points interpolated from simulated FID, up to maximum simulated timepoint (Tmax)
            List<Point> interpolatedTransferFunction = new List<Point>();

            // array of interpolated points, which includes the simulated FID (up to Tmax)
            // plus additional zero filled points to cover the entire requested time domain points
            Point[] interpolatedTransferFunctionWithZeroFilling = new Point[bS.TimeDomain];
            
            // iterate through each desired timepoint,
            // creating a new point at the interpolated amplitude
            // or zero filling for each point in the FID
            for (int k = 0; k < bS.TimeDomain; k++)
            {
                // the current time value is calculated based on the time increment (1/Spectral Width)
                double timepoint = (k + 1) * bS.Increment;

                // Interpolation will only work if the calculated timepoint is within the maximum timepoint
                // sampled from the Spinach data
                if (timepoint <= simData.Tmax)
                {
                    interpolatedTransferFunction.Add(new Point(timepoint, realDataInterpolation.Interpolate(timepoint), imagDataInterpolation.Interpolate(timepoint)));
                }
                else
                {
                    // all remaining points in time domain are larger than sampled points from spinach
                    // (i.e. larger than Tmax), thus we want to zero fill for the rest of the TimeDomain
                    // points, then exit interpolation.

                    interpolatedTransferFunction.CopyTo(interpolatedTransferFunctionWithZeroFilling);

                    for (int i = k; i < bS.TimeDomain; i++)
                    {
                        timepoint = (i + 1) * bS.Increment;
                        interpolatedTransferFunctionWithZeroFilling[i] = new Point(timepoint, 0, 0);
                    }

                    break;
                }
            }

            // ALISTAIR 06/03/2012: Add interpolated data to list for further interpolation across tauC.
            // Exclude simulated FIDs for very small values of TauC (TauC < 1 ns), since these do not reflect
            // the relaxation properties that are being modelled in the interpolation
            if (simData.TauC >= 1)
            {
                tauCsForInterpolation.Add(simData.TauC);
                timeDomainInterpolatedData.Add(interpolatedTransferFunction);
            }

            // Diagnostic output for debugging
            /*
            String diagnosticOutput = "Interpolated Transfer Function:\n";
            diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
            diagnosticOutput += "Datapoint,Time,Value\n";
            for (int k = 0; k < InterpolatedTransferFunction.Count; k++)
            {
                diagnosticOutput += (k + 1) + "," + InterpolatedTransferFunction.ElementAt(k).T.ToString()
                                    + "," + InterpolatedTransferFunction.ElementAt(k).V.Real.ToString() + "\n";
            }
            System.Diagnostics.Debug.WriteLine(diagnosticOutput);
            */

            // ALISTAIR 12/03/2012: For interpolation where no zero filling is required (Tmax > sampling time)
            // then directly copy interpolated function data to array of zero filling data (no zero filling required)
            if (interpolatedTransferFunction.Count() == bS.TimeDomain)
            {
                interpolatedTransferFunction.CopyTo(interpolatedTransferFunctionWithZeroFilling);
            }

            // Diagnostic output for debugging
            /*
            String diagnosticOutput = "Interpolated Transfer Function With Zero Filling:\n";
            diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
            diagnosticOutput += "Datapoint,Time,Value\n";
            for (int k = 0; k < interpolatedTransferFunctionWithZeroFilling.Length; k++)
            {
                diagnosticOutput += (k + 1) + "," + interpolatedTransferFunctionWithZeroFilling.ElementAt(k).T.ToString()
                                    + "," + interpolatedTransferFunctionWithZeroFilling.ElementAt(k).V.Real.ToString() + "\n";
            }
            System.Diagnostics.Debug.WriteLine(diagnosticOutput);
            */

            // Apply apodization function over number of points specified in the time domain
            // Using cosine bell apodization function (ssb=2, SIN or QSIN)
            double[] apodizationFunction = new double[bS.TimeDomain];
            for (int k = 0; k < bS.TimeDomain; k++)
            {
                apodizationFunction[k] = Math.Cos((k + 1) * ((Math.PI / 2) / bS.TimeDomain));
            }
            
            // Diagnostic output for debugging
            /*
            String diagnosticOutput = "Apodization Function:\n";
            diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
            diagnosticOutput += "Timepoint,Value\n";
            for (int k = 0; k < apodizationFunction.Length; k++)
            {
                diagnosticOutput += (k + 1) + "," + apodizationFunction.ElementAt(k).ToString() + "\n";
            }
            System.Diagnostics.Debug.WriteLine(diagnosticOutput);
            */

            // ALISTAIR: Comment added 19/03/2012
            // Since any individually calculated signal peak will also contain noise, a number of of trials are
            // required to average out the noise (which will average to zero) from the signal
            // TODO: The number of trials required to get a reliable signal will depend on the initial S-to-N
            // from the 1D experiment. If the percentage of noise is high (low signal-to-noise), then lots of
            // trials will be required to get a reliable signal; if the percentage of noise is low
            // (high signal-to-noise), then fewer trial will be required to get a reliable signal.
            // Need to define a threshold at which the signal is considered to be within an acceptable
            // margin of error (e.g +/- 1% error)
            // int numOfTrials = 1/bS.OneDimSigToNoise * SomeFunction
            //int numOfTrials = 64 * 64;
            int numOfTrials = 8;           
            Random RandomGenerator = new Random();

            // ALISTAIR: COMMENT ADDED 24/01/2011
            // Generate an array of 'trials' numbers for each
            // parameter of interest
            double[] PeakHeight = new double[numOfTrials];
            double[] NoiseRMSD = new double[numOfTrials];
            double[] SN = new double[numOfTrials];
            double[] Width = new double[numOfTrials];

            // ALISTAIR: COMMENT ADDED 24/01/2011
            // set the number of points in frequency domain
            // instead use bS.NoOfPointsSI;
            //int SI = 2048;
            //int SI = 8192;
            //int SI = 512;

            Complex[] FIDwithSIpoints = null;
            Complex[] spectrumPoints = null;

            // ALISTAIR: COMMENT ADDED 24/01/2011
            // Iterate through all 'trials'
            for (int m = 0; m < numOfTrials && !this.Terminated; m++)
            {
                // Random Noise Function over number of point in time domain
                // ALISTAIR Modified 18/08/2011
                // Rather than using arbitrary 10% noise (0.1),
                // generate noise inversely proportional to signal-to-noise ratio from 1D experiment
                // (parameter input by user)
                // The noise applied is the inverse of signal-to-noise from 1D (for all transients)
                // (e.g. a 1D S-to-N of 10 applies a noise of 1/10 = 0.1
                //       a 1D S-to-N of 20 applies a noise of 1/20 = 0.05)
                // Since the total noise is spread over both the real and imaginary signals,
                // the noise added to each is divided by square root of 2
                double noiseToSignal = 1 / bS.OneDimSigToNoise;
                //S/N (over all transients) should be divided the square root of the
                //number of transients (scans) to get the S/N per scan, which can be applied in the
                //time domain for the indirect dimension
                //double OneDimSignalToNoisePerScan = bS.OneDimSigToNoise / Math.Sqrt(bS.NoOfTransients);
                //double noiseToSignal = 1 / OneDimSignalToNoisePerScan;
                Complex32[] RandomNoiseFunction = new Complex32[bS.TimeDomain];

                for (int k = 0; k < bS.TimeDomain; k++)
                {
                    // ALISTAIR 14/03/2012: 
                    //RandomNoiseFunction[k] = new Complex32((float)RandomGenerator.NextDouble(), (float)noiseToSignal * (float)RandomGenerator.NextDouble());
                    // Need to acquire separate random numbers for the real and imaginary components
                    float randomNoiseReal = (float)((noiseToSignal * (RandomGenerator.NextDouble() - 0.5)) / Math.Sqrt(2)); 
                    float randomNoiseImag = (float)((noiseToSignal * (RandomGenerator.NextDouble() - 0.5)) / Math.Sqrt(2)); 
                    RandomNoiseFunction[k] = new Complex32(randomNoiseReal, randomNoiseImag);
                }

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // For each trial, FID is an array of complex numbers, the length being the number of points
                // required after Fourier transform into the frequency domain
                // Note that since SI (points in freq domain) is likely to be larger than TD (points in time domain),
                // so further zero filling will take place for all points larger than TD
                // Modified: Removed R2Effective from calculation
                // Modified: To cope with both real and imag components
                FIDwithSIpoints = new Complex[bS.NoOfPointsSI];
                for (int k = 0; k < bS.TimeDomain && k < bS.NoOfPointsSI; k++)
                {
                    //FID[k] = (InterpolatedTransferFunction[k].V * Math.Exp(-bS.R2Effective * InterpolatedTransferFunction[k].T) + RandomNoiseFunction[k]) * apodizationFunction[k];
                    //FIDwithSIpoints[k] = (interpolatedTransferFunctionWithZeroFilling[k].V.Real + RandomNoiseFunction[k].Real) * apodizationFunction[k];
                    FIDwithSIpoints[k].Real = (interpolatedTransferFunctionWithZeroFilling[k].V.Real + RandomNoiseFunction[k].Real) * apodizationFunction[k];
                    FIDwithSIpoints[k].Imag = (interpolatedTransferFunctionWithZeroFilling[k].V.Imaginary + RandomNoiseFunction[k].Imaginary) * apodizationFunction[k];
                }

                // Diagnostic output for debugging
                /*
                // Only sample first two iterations of 'trials'
                if (m < 2)
                {
                    String diagnosticOutput = "InterpolatedFIDWithNoiseAndApodization\n";
                    diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
                    diagnosticOutput += "Trial: " + (m+1) + "\n";
                    diagnosticOutput += "Timepoint,Real Value,Imag Value\n";
                    for (int k = 0; k < FIDwithSIpoints.Length; k++)
                    {
                        diagnosticOutput += (k + 1) + "," + FIDwithSIpoints.ElementAt(k).Real.ToString()
                                         + "," + FIDwithSIpoints.ElementAt(k).Imag.ToString() + "\n";
                    }
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */


                // ALISTAIR 02/03/2012: perform a forward Fourier transform on the FID data
                // The input array is modified by the TransformForward function
                ComplexFourierTransformation cft = new ComplexFourierTransformation(TransformationConvention.Matlab);
                cft.TransformForward(FIDwithSIpoints);
                // copy the Fourier transformed FID points as an array of complex numbers representing a spectrum 
                spectrumPoints = FIDwithSIpoints;
                 
                // Diagnostic output for debugging
                /*
                // Only sample first two iterations of 'trials' for TauC = 10 ns
                if (m < 2 && simData.TauC == 10)
                {
                    String diagnosticOutput = "AfterForwardFourierTransformation\n";
                    diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
                    diagnosticOutput += "Trial: " + (m+1) + "\n";
                    diagnosticOutput += "Timepoint,Real Value,Imag Value\n";

                    for (int k = 0; k < FID.Length; k++)
                    {
                        diagnosticOutput += (k + 1) + "," + FID.ElementAt(k).Real.ToString() + "," + FID.ElementAt(k).Imag.ToString() + "\n";
                    }

                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */

                // FFT Shift - change first and second halfs of the spectrum
                int middle = bS.NoOfPointsSI / 2;
                for (int j = 0, k = middle; j < middle && k < bS.NoOfPointsSI; j++, k++)
                {
                    Complex temp = spectrumPoints[j];
                    spectrumPoints[j] = spectrumPoints[k];
                    spectrumPoints[k] = temp;
                }

                /***************************************************
                 * 
                 * THE PLAN FOR CREATING SPECTRUM WITH FEQUENCIES SHIFTED FOLLOWING A BELL-SHAPED DISTRIBUTION OF DENSITIES
                 * (ALL WITH THE SAME AMPLITUDE)
                 * 
                 * Aim. Should have a total of N peaks, where n is determined by number of peaks expected, depending on the sampling scheme.
                 *
                 * 1. Display diagnostics below (after switching two halves of spectrum), to output data from FT process.
                 *    Analyse used excel.
                 * 2. Find appropriate reverse FT function and apply this to the FTed data.
                 * 3. FT this again, and analyse using excel again, to make sure that the FT and reverse FT processes result in the same data.
                 * 4. Need to freq. shift the FTed data (preferably using a circular freq. shift method), perhaps using a for loop to shift
                 *    each data point along a few points in the FTed array, then moving the last few points in the array to the beginning of the array
                 *    (perhaps create another array - or a 2D array, with each dataset being freq shifted from the previous in the array).
                 * 5. Display diagnostics of the new array to make sure it has visibly fequency shifted.
                 * 6. Add the values of the old array and the new array together, to create a third array that cadjuontains two peaks.
                 * 7. Display diagnostics to confirm two frequencies are displayed.
                 * 8. Repeat the process in a simple loop to create about 10 frequencies, initially just shifted equally from each other.
                 * 9. Display diagnostics to confirm all frequencies are displayed.
                 * 10. Perhaps adjust so that original freq is in centre, at 0 Hz.
                 * 11. Need to adjust to N, where N is the expected number of peaks.
                 * 12. Need to reverse FT each peak separately then add values together (or reverse FT of spectrum with all peaks combined)??
                 * 13. Final interferogram will be the addition of all contributions from each peak in the Fourier transformed spectrum
                 *      (but reverse Fourier transformed so that the data is in the time domain)
                 * 14. Send data to ALNMR - need to work out required format.
                 * 
                 * NEXT, CALCULATE THE POINTS NECESSARY FOR NON-LINEAR SAMPLING
                 * 
                 * 15. Get all the points (in time domain?) for the 15N dimension. Plot using excel.
                 * 16. Get all the points (in time domain?) for the 13C dimension. Plot using excel.
                 * 17. Extend to two dimensional array of points in 15N, 13C dimensions.
                 *     (using a simple mathematical operation - just multiply values together???)
                 * 18. Plot and see if this makes any sense!!
                 * 19. Start implementing some sort of sampling scheme, e.g. take multiple sample in different directions
                 *     from origin.
                 * 20. Send non-linear sample data to ALNMR, see what it spits out. Any way to verify the results?
                 * 
                 * *************************************************/

                // Diagnostic output for debugging
                /*
                // Only sample first two iterations of 'trials' for TauC = 10 ns
                if (m < 2 && simData.TauC == 10)
                {
                    String diagnosticOutput = "AfterSwitchingTwoHalvesOfSpectra\n";
                    diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
                    diagnosticOutput += "Trial: " + (m + 1) + "\n";
                    diagnosticOutput += "Timepoint,Real Value,Imag Value\n";

                    for (int k = 0; k < spectrumPoints.Length; k++)
                    {
                        diagnosticOutput += (k + 1) + "," + spectrumPoints.ElementAt(k).Real.ToString() + "," + spectrumPoints.ElementAt(k).Imag.ToString() + "\n";
                    }

                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */

                // BaseLine is a arithmetic mean of first quarter points
                // ALISTAIR: COMMENT ADDED 24/01/2011
                // This is making an assumption that the desired
                // peak will not be within the first quarter of the spectrum
                // i.e. the first quarter is all baseline signal (noise)
                Complex baseLine = new Complex(0, 0);

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // each data point in the first quarter of the spectrum
                // is added to accumulate the baseline data
                int quarterSI = bS.NoOfPointsSI / 4;
                for (int p = 0; p < quarterSI; p++)
                {
                    baseLine += spectrumPoints[p];
                }

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // baseline becomes an average of all of the points in the first quarter
                baseLine /= quarterSI;

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // Substract the baseline from each point in the spectrum
                // i.e. baseline correction
                for (int k = 0; k < bS.NoOfPointsSI; k++)
                {
                    spectrumPoints[k] -= baseLine;
                }

                // Diagnostic output for debugging
                /*
                // Only sample first two iterations of 'trials' for TauC = 10 ns
                if (m < 2 && simData.TauC == 10)
                {
                    String diagnosticOutput = "AfterBaselineCorrection\n";
                    diagnosticOutput += "TauC: " + simData.TauC.ToString() + "\n";
                    diagnosticOutput += "Trial: " + (m + 1) + "\n";
                    diagnosticOutput += "Timepoint,Real Value,Imag Value\n";

                    for (int k = 0; k < spectrumPoints.Length; k++)
                    {
                        diagnosticOutput += (k + 1) + "," + spectrumPoints.ElementAt(k).Real.ToString() + "," + spectrumPoints.ElementAt(k).Imag.ToString() + "\n";
                    }

                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // For each trial, the peak height becomes the max real vaule in the spectrumPoints array
                PeakHeight[m] = spectrumPoints.Max(ff => ff.Real);
                //  Complex baseLine = spectrumPoints.Take(SI / 4).Sum(pp => pp) / SI / 4; // Sum does not work for complex

                // Diagnostic output for debugging
                /*
                // Only sample first two iterations of TauC (i.e. <= 10 ns)
                if (simData.TauC <= 10)
                {
                    String diagnosticOutput = "TauC,Trial,Peak Height (Real),";
                    diagnosticOutput += simData.TauC.ToString() + "," + (m + 1) + "," + PeakHeight[m];
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // For each trial, calculate noise as perturbation of the baseline.
                // (since the expected average of noise is zero, the magnitude
                //  represents the deviation from this expected value)
                // 1. For each point of the baseline, multiply the complex number by its
                //    complex conjugate, so that a positive real number is produced (squared deviation, SD)
                // 2. Sum all the squared noise values (sum of squared deviation)
                // 3. Divide the sum by the number of baseline points sampled (mean of squared deviation, MSD)
                // 4. The RMSD of the noise for each trial becomes the square root of the mean value
                double sumOfSquares = 0;
                double meanOfSquares = 0;
                for (int p = 0; p < quarterSI; p++)
                {
                    // get the sum of the squares
                    // complex * complex conjugate is effectively the square of a complex number
                    sumOfSquares += (spectrumPoints[p] * spectrumPoints[p].Conjugate).Real; // c*c.Conjugate is real
                    // sumOfSquares += spectrumPoints[p].Square().Real;
                }
                meanOfSquares = sumOfSquares / quarterSI;
                NoiseRMSD[m] = Math.Sqrt(meanOfSquares);


                // Diagnostic output for debugging
                /*                
                // Only sample first two iterations of TauC (i.e. <= 10 ns)
                if (simData.TauC <= 10)
                {
                    String diagnosticOutput = "TauC,Trial,NoiseRMSD,";
                    diagnosticOutput += simData.TauC.ToString() + "," + (m + 1) + "," + NoiseRMSD[m];
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // The peak signal is the max of the spectrum values (real components)
                // divided by current noise to get Signal-to-Noise value
                SN[m] = PeakHeight[m] / NoiseRMSD[m];

                // Diagnostic output for debugging
                /*
                // Only sample first two iterations of TauC (i.e. <= 10 ns)
                if (simData.TauC <= 10)
                {
                    String diagnosticOutput = "TauC,Trial,Signal-To-Noise,";
                    diagnosticOutput += simData.TauC.ToString() + "," + (m + 1) + "," + SN[m];
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */

                // ALISTAIR: Modified 16/03/2012
                // To calculate the linewidth of the peak:
                // 1. Assume peak is centred in middle of spectrum (which it usually is from Spinach simulations)
                // 2. Find the "real" value associated with peak (PeakHeight[m])
                // 3. Find the "half height" of the peak, which will be PeakHeight[m]/2
                // 4. Traverse through the points of the spectrum, from the middle downwards, and find the
                //    first "real" value (spectrumPoints[k].Real) that is less than or equal to the half height (PeakHeight[m]/2).
                //    Once found, this value of k represents a point in the spectrum at which the peak is at half height
                //    (the point at half height below the peak maximum)
                // 5. The middle point of the spectrum (SI/2) represents the peak point
                // 6. Half width of half height (in number of points) is given by midpoint - point at half height (SI/2 - k)
                // 7. There will be a corresponding half height point on the upper part of the spectrum, so need
                //    to multiply by two to get full width (in number of points) at half height (2*(SI/2 - k))
                // 8. To convert linewidth to a frequency, first find the spectral width (SW, in Hz) = 1/ bS.increment
                // 9. The width per point (in Hz) = SW / number of points (SI), i.e. (1/bS.Increment)/SI
                // 10. The linewidth (in Hz) = full width at half height (in number of points) * width per point (in Hz),
                //     i.e. (2*(SI/2 - k)) * ((1/bS.Increment)/SI)
                for (int k = middle; k >= 0; k--)
                {
                    if (spectrumPoints[k].Real <= (PeakHeight[m] / 2))
                    {
                        Width[m] = (2 * ((bS.NoOfPointsSI / 2) - k)) * ((1 / bS.Increment) / bS.NoOfPointsSI);
                        break;
                    }
                }

                // Diagnostic output for debugging
                /*
                // Only sample first two iterations of TauC (i.e. <= 10 ns)
                if (simData.TauC <= 10)
                {
                    String diagnosticOutput = "TauC,Trial,Width,";
                    diagnosticOutput += simData.TauC.ToString() + "," + (m + 1) + "," + Width[m];
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                }
                */

            }  // end of trials loop

            if (this.Terminated) return null;

            // ALISTAIR: COMMENT ADDED 24/01/2011
            // Create a new BenchmarkResult object for displaying output
            BenchmarkResult bR = new BenchmarkResult();

            // ALISTAIR: COMMENT ADDED 24/01/2011
            // Copy Dimension, Time Domain and TauC parameters directly from BenchmarkSetup object
            bR.Dimension = bS.Dimension;
            bR.TimeDomain = bS.TimeDomain;
            bR.TauC = simData.TauC;

            // ALISTAIR: MODIFIED 27/03/2012
            // Averaged Peak Height and RMSD
            // To calculate RMSD of peak height (and all other RMSD below):
            // 1. subtract each individual peak height value from the average (i.e deviation from mean).
            // 2. Only interested in magnitude of the difference from the mean,
            //    so square the value (squared deviation, SD)
            // 3. Get the sum of all differences from the mean (sum of squared deviation)
            // 4. Divide by the number of 'trials' (mean of squared deviation, MSD)
            // 5. The RMSD becomes the square root of the mean
            double signal_avg = PeakHeight.Sum() / numOfTrials;
            double signal_rmsd = Math.Sqrt(PeakHeight.Sum(ph => Math.Pow(signal_avg - ph, 2)) / numOfTrials);
            bR.Signal = new PhysicalValue<double>(signal_avg, signal_rmsd);

            // ALISTAIR: MODIFIED 27/03/2012
            // Averaged NoiseRMSD and the "RMSD of the noise RMSD"
            // i.e. for each trial, the RMSD of noise is calculated for the individual spectrum;
            // over all trials, the RMSD is calculated across all spectra, with respect to the
            // individually calculated noise RMSD for each spectrum.
            double noise_avg = NoiseRMSD.Sum() / numOfTrials;
            double noise_rmsd = Math.Sqrt(NoiseRMSD.Sum(noise => Math.Pow(noise_avg - noise, 2)) / numOfTrials);
            bR.Noise = new PhysicalValue<double>(noise_avg, noise_rmsd);

            // ALISTAIR: MODIFIED 27/03/2012
            // Averaged Linewidth and RMSD
            // The linewidth of the peak can never be more precise than the frequency between adjacent
            // points on the discretized spectrum
            // If this is the case, re-assign the error to be the frequency between points in the spectrum
            double width_avg = Width.Sum() / numOfTrials;
            double width_rmsd = Math.Sqrt(Width.Sum(wi => Math.Pow(width_avg - wi, 2)) / numOfTrials);
            if (width_rmsd < (1 / bS.Increment) / bS.NoOfPointsSI)
            {
                width_rmsd = (1 / bS.Increment) / bS.NoOfPointsSI;
            }

            bR.Linewidth = new PhysicalValue<double>(width_avg, width_rmsd);

            // ALISTAIR: MODIFIED 27/03/2012
            // Averaged S/N and RMSD
            double sn_avg = SN.Sum() / numOfTrials;
            double sn_rmsd = Math.Sqrt(SN.Sum(sn => Math.Pow(sn_avg - sn, 2)) / numOfTrials);
            bR.SignalToNoise = new PhysicalValue<double>(sn_avg, sn_rmsd);

            bR.StNtTD = sn_avg / bS.TimeDomain;

            // ALISTAIR 27/03/2012:
            // Temporarily add additional values, for testing purposes
            /*
            bR.StoNfrom1D = bS.OneDimSigToNoise;
            bR.SI = bS.NoOfPointsSI;
            bR.NoOfScans = bS.NoOfTransients;
            bR.StoNdivRootTD = sn_avg / Math.Sqrt(bS.TimeDomain);
            */

            // Diagnostic output for debugging
            /*
            String diagnosticOutput = "Calculated statistics\n\n";
            diagnosticOutput += "Dimension," + bR.Dimension + "\n";
            diagnosticOutput += "TimeDomain," + bR.TimeDomain + "\n";
            diagnosticOutput += "TauC," + simData.TauC.ToString() + "\n\n";
            diagnosticOutput += "Average Peak Height," + signal_avg + "\n";
            diagnosticOutput += "Peak Height RMSD," + signal_rmsd + "\n";
            diagnosticOutput += "Max Signal (bR.Signal) Value," + bR.Signal.Value + "\n";
            diagnosticOutput += "Max Signal (bR.Signal) Error (RMSD)," + bR.Signal.Error + "\n\n";
            diagnosticOutput += "Average Noise," + noise_avg + "\n";
            diagnosticOutput += "Noise RMSD," + noise_rmsd + "\n";
            diagnosticOutput += "Noise (bR.Noise) Value," + bR.Noise.Value + "\n";
            diagnosticOutput += "Noise (bR.Noise) Error (RMSD)," + bR.Noise.Error + "\n\n";
            diagnosticOutput += "Average Linewidth," + width_avg + "\n";
            diagnosticOutput += "Linewidth RMSD," + width_rmsd + "\n";
            diagnosticOutput += "Linewidth (bR.Linewidth) Value," + bR.Linewidth.Value + "\n";
            diagnosticOutput += "Linewidth (bR.Linewidth) Error (RMSD)," + bR.Linewidth.Error + "\n\n";
            diagnosticOutput += "Average S-to-N," + sn_avg + "\n";
            diagnosticOutput += "S-to-N RMSD," + sn_rmsd + "\n";
            diagnosticOutput += "S-to-N (bR.SignalToNoise) Value," + bR.SignalToNoise.Value + "\n";
            diagnosticOutput += "S-to-N (bR.SignalToNoise) Error (RMSD)," + bR.SignalToNoise.Error + "\n\n";
            diagnosticOutput += "S-to-N per TD (bR.StNtTD)," + bR.StNtTD + "\n";
            System.Diagnostics.Debug.WriteLine(diagnosticOutput);
            */

            // ALISTAIR: COMMENT ADDED 24/01/2011
            // Construct the 1D spectrum to be displayed as an array of points
            Point[] spectrum = new Point[bS.NoOfPointsSI];
            for (int i = 0; i < bS.NoOfPointsSI; i++)
            {
                // ALISTAIR: COMMENT ADDED 23/03/2012
                // For each new Point:
                // note that since 'T' is now in the frequency domain, this represents frequency rather than time.
                // the frequency of the point is calculated by:
                // 1. spectral width (SW, Hz) = 1 / bS.Increment
                // 2. width per point (Hz)  = SW / number of points = (1 / bS.Increment) / bS.NoOfPointsSI
                // 3. freq of point (Hz) = point * width per point = (i+1) * (1 / bS.Increment) / bS.NoOfPointsSI
                //                                                 = (i+1) / bS.Increment / bS.NoOfPointsSI
                // 4. since we want to calibrate the spectrum so that the zero point (0 Hz) is at the centre of the peak,
                //    which is in the centre of the spectrum, rather than the lowest point of the spectrum, we must subtract
                //    half of the width of the spectrum from each point = SW / 2 = (1/bS.Increment) * 0.5 = 0.5 / bS.Increment
                // 5. calibrated freq of point (Hz) = (i+1) / bS.Increment / bS.NoOfPointsSI - 0.5 / bS.Increment
                // The value, 'V', becomes a combination of the Real and Imag portions of the spectrum value.
                Point p = new Point((i + 1) / bS.Increment / bS.NoOfPointsSI - 0.5 / bS.Increment, spectrumPoints[i].Real, spectrumPoints[i].Imag);
                spectrum[i] = p;
            }

            bR.Spectrum = spectrum;

            return bR;

        }  // End CalculateBenchmarks()

        // ALISTAIR: COMMENT ADDED 24/01/2011
        // Executes analysis of the benchmark data for all dimensions
        private void RunEverything()
        {
            DateTime starting = DateTime.Now;
            this.ExecutionStatus = "Begin Execution";

            try
            {
                // flag indicating current dimension in list of dimensions (InputData)
                int dimensionNumber = 0;

                // ALISTAIR: COMMENT ADDED 24/01/2011
                // Traverse through each <BenchmarkSetup> set of data
                // i.e. each dimension, calculating statistics
                foreach (var bS in this.InputData)
                {
                    // increment dimension number
                    dimensionNumber++;

                    // ALISTAIR 13/03/2012: initialise the interpolation lists for each indirect dimension
                    tauCsForInterpolation = new List<double>();
                    timeDomainInterpolatedData = new List<List<Point>>();

                    DateTime beg = DateTime.Now;
                    this.ExecutionStatus = string.Format("Processing Dimension {0}.", bS.Dimension);

                    // ALISTAIR 16/08/2011:
                    // Updated to perform statistics for each set of 'Real' data in each dimension
                    foreach (var simData in bS.DimensionData)
                    {
                        BenchmarkResult bM = this.CalculateBenchmarks(bS, simData);
                        if (bM == null) break;

                        // ALISTAIR: COMMENT ADDED 24/01/2011
                        // Record the execution time for each set of results
                        bM.ExecutionTime = DateTime.Now - beg;
                        this.Results.Add(bM);
                    }

                    if (this.Terminated) break;

                    // ALISTAIR Added 23/05/2011
                    // In order to provide analytic data for any desired rotational correlation time (tauC),
                    // it is necessary to interpolate the simulated tauC data.
                    // First: Interpolate from sparse sampling of Spinach to the number of points specified by the Time Domain parameter
                    //        (already performed in CalculateBenchmarks() method as part of the analytics calculations)
                    // Second: For array of all Time Domain points, interpolate in the "tauC dimension" for the user specified tauC value.
                    // This 2D interpolation is calculated independently for each indirect dimension.
                    // For example:
                    // In Spinach: Sampled 50 timepoints for tauC = 10 ns to tauC = 60 ns, in steps of 10 ns.
                    // In NMRplus: Time Domains = 200 is supplied as a parameter.
                    // 1st Dimension: For each of the 6 provided values of TauC, interpolate to extend from the 50 timepoints to 200 timepoints.
                    // 2nd Dimension: For each of the 200 timepoints, interpolate using the known tauC values for a the user's tauC

                    // Declare a 2D list of points which will transpose the time domain interpolated data,
                    // (which is a list of sampled tauC values each containing a list of timepoint/value pairs),
                    // into a list of timepoints each containing a list of points (one for each tauC).
                    // This is used to facilitate interpolation across tauC 
                    List<List<Point>> tauCDataBefore2ndInterpolation = new List<List<Point>>();

                    // Diagnostic output for debugging
                    /*
                    String diagnosticOutput = "List of Time Domain Interpolated TauC Values\n";
                    diagnosticOutput += "timeDomainInterpolatedData.Count = " + timeDomainInterpolatedData.Count + "\n";
                    foreach (List<Point> pointsForTimeDomainTuaC in timeDomainInterpolatedData)
                    {
                        diagnosticOutput += "pointArray.Length," + pointsForTimeDomainTuaC.Count;
                        diagnosticOutput += "\nTimepoint,Real Value,Imag Value\n";
                        foreach (Point aTimeInterpolatedPoint in pointsForTimeDomainTuaC)
                        {
                            diagnosticOutput += aTimeInterpolatedPoint.T.ToString() + "," +
                                                              aTimeInterpolatedPoint.V.Real.ToString() + "," +
                                                              aTimeInterpolatedPoint.V.Imaginary.ToString() + "\n";
                        }
                    }
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                    */

                    // the initial amplitude of the FID at time zero,
                    // which will be different for each tauC
                    // ALISTAIR Modified 18/10/2012:
                    // Changed initial amplitude to an array the size of the number of simulated tauCs
                    double[] initialAmpReal = new double[tauCsForInterpolation.Count];
                    double[] initialAmpImag = new double[tauCsForInterpolation.Count];
                    // (an initial amplitude of 1 can be assumed whenever we wish to observe
                    // relaxation effects only, independent of loss of magnetisation during pulse sequence)
                    //double initialAmpReal = 1;

                    // transpose the 2D list of FIDs interpolated for the desired timepoints
                    // into a 2D list prepared for interpolation across tauC
                    for (int j=0; j < timeDomainInterpolatedData.Count; j++)
                    {
                        // Due to loss of signal during magnetisation transfer, the initial amplitude will be less than one,
                        // and the magnitude of the loss of signal will be proportional to the rotational correlation time (tauC).
                        // The initial amp can be taken as the amplitude at the first timepoint of each tauC simulation
                        initialAmpReal[j] = timeDomainInterpolatedData[j][0].V.Real;
                        initialAmpImag[j] = timeDomainInterpolatedData[j][0].V.Imaginary;

                        for (int k = 0; k < timeDomainInterpolatedData[j].Count; k++)
                        {
                            if (j == 0)
                            {
                                tauCDataBefore2ndInterpolation.Add(new List<Point>());
                            }

                            tauCDataBefore2ndInterpolation[k].Add(timeDomainInterpolatedData[j][k]);
                        }
                    }

                    // ALISTAIR Added 18/10/2012:
                    // Create an interpolation function to calculate the initial amplitude for user input tauC value
                    // Using natural cubic spline interpolation
                    // The first parameter is the list of simulated tauC values
                    // The second parameter is the list of initial values associated with each simultaed tauC
                    // Note: Since the interpolation method can not handle complex numbers,
                    // the real and imaginary components of each complex number are interpolated independently
                    IInterpolationMethod realInitialAmpInterpolation = Interpolation.CreateNaturalCubicSpline(tauCsForInterpolation, initialAmpReal);
                    IInterpolationMethod imagInitialAmpInterpolation = Interpolation.CreateNaturalCubicSpline(tauCsForInterpolation,initialAmpImag);
                    double initialRealAmpForUserTauC = realInitialAmpInterpolation.Interpolate(userTauC);
                    double initialImagAmpForUserTauC = imagInitialAmpInterpolation.Interpolate(userTauC);

                    // Diagnostic output for debugging
                    /*                    
                    String diagnosticOutput = "TauC Data Before Interpolation (Transpose of Time Domain Interpolated)\n";
                    diagnosticOutput += "tauCDataBefore2ndInterpolation.Count = " + tauCDataBefore2ndInterpolation.Count + "\n";
                    for (int i=0; i<tauCsForInterpolation.Count; i++)
                    {
                        diagnosticOutput += "initialAmpReal[" + i + "] = " + initialAmpReal[i] + "\n";
                    }
                    foreach (List<Point> pointsInTauCBeforeInterpolation in tauCDataBefore2ndInterpolation)
                    {
                        diagnosticOutput += "pointArray.Length," + pointsInTauCBeforeInterpolation.Count;
                        diagnosticOutput += "\nTimepoint,Real Value,Imag Value\n";
                        int s = 0;
                        foreach (Point aTimeInterpolatedPoint in pointsInTauCBeforeInterpolation)
                        {
                            diagnosticOutput += aTimeInterpolatedPoint.T.ToString() + "," +
                                                              aTimeInterpolatedPoint.V.Real.ToString() + "," +
                                                              aTimeInterpolatedPoint.V.Imaginary.ToString();
                            diagnosticOutput += ",TauC: " + tauCsForInterpolation.ElementAt(s) + "\n";
                            s++;
                        }
                    }                
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                    */                  

                    // array of points representing a simulated FID at the user specified tauC
                    Point[] userTauCInterpolatedPoints = new Point[bS.TimeDomain];

                    // Loop through each timepoint, using the known points from simulated tauCs
                    // to interpolate a point for the user tauC
                    // new points are calculated based on the formula:
                    // amplitude(t) = amplitude(0) * e-(C*tauC*t)
                    // where C is a constant that can be calculated for each timepoint based on the
                    // simulated points of each tauC
                    // interpolation of real and imaginary components are performed independently
                    // Programmatically, this can be written as:
                    // amplitudeReal = initialAmpReal * exp^(-1*constantCReal*currentTauC*currentTimepoint)
                    // amplitudeImag = initialAmpImag * exp^(-1*constantCImag*currentTauC*currentTimepoint)
                    for (int timepointIndex = 0; timepointIndex < tauCDataBefore2ndInterpolation.Count; timepointIndex++)
                    {
                        // get the current timepoint, which is the same for all tauC data interpreted in the time domain,
                        // so can extract the value from the first simulated tauC
                        double currentTimepoint = tauCDataBefore2ndInterpolation[timepointIndex].First().T;

                        // initialize variables for the sum and average of a constant
                        // required for tauC interpolation formula
                        double sumOfCForSimulatedTauCsReal = 0;
                        double avgConstantCReal = 0;
                        double sumOfCForSimulatedTauCsImag = 0;
                        double avgConstantCImag = 0;

                        // initialize the amplitude for the real and imag components at the user tauC
                        double amplitudeReal = 0;
                        double amplitudeImag = 0;

                        // Iterate through all tauC for which a simulated datapoint is available at the current timepoint.
                        // Based on the knonwn datapoints, an average value for constant C can be calculated.
                        for (int tauCIterator = 0; tauCIterator < tauCDataBefore2ndInterpolation[timepointIndex].Count; tauCIterator++)
                        {
                            // get the current tauC value and associated datapoint
                            double currentTauC = tauCsForInterpolation[tauCIterator];
                            double realValueForTauC = tauCDataBefore2ndInterpolation[timepointIndex][tauCIterator].V.Real;
                            double imagValueForTauC = tauCDataBefore2ndInterpolation[timepointIndex][tauCIterator].V.Imaginary;

                            // Apply current tauC, current time, and current real value (amplitude) to get a value for the constant C.
                            // Rearrange formula amplitudeReal = initialAmp * exp^(-1*constantC*currentTauC*currentTimepoint)
                            // exp^(-1*constantC*currentTauC*currentTimepoint) = amplitudeReal/initialAmp
                            // -1*constantC*currentTauC*currentTimepoint = logBaseE(amplitudeReal/initialAmp)
                            // constantC*currentTauC*currentTimepoint = -1*logBaseE(amplitudeReal/initialAmp)
                            // constantC = (-1*logBaseE(amplitudeReal/initialAmp))/(currentTauC*currentTimepoint)
                            double constantCForCurrentTauCReal = (-1 * Math.Log(realValueForTauC / initialAmpReal[tauCIterator])) / (currentTauC * currentTimepoint);
                            // ALISTAIR Added 18/10/2012:
                            // Add check for imag data = zero before calculation
                            double constantCForCurrentTauCImag = 0;
                            if (initialAmpImag[tauCIterator] != 0)
                            {
                                constantCForCurrentTauCImag = (-1 * Math.Log(imagValueForTauC / initialAmpImag[tauCIterator])) / (currentTauC * currentTimepoint);
                            }

                            // add the calculated constant C to a sum for all simulated tauCs
                            sumOfCForSimulatedTauCsReal += constantCForCurrentTauCReal;
                            sumOfCForSimulatedTauCsImag += constantCForCurrentTauCImag;
                        }

                        // calculate the average constant C based on calculations from all tauC datapoints for this timepoint
                        // (the more datapoints that are available, the more accurate the constant C will become)
                        avgConstantCReal = sumOfCForSimulatedTauCsReal / tauCDataBefore2ndInterpolation[timepointIndex].Count;
                        avgConstantCImag = sumOfCForSimulatedTauCsImag / tauCDataBefore2ndInterpolation[timepointIndex].Count;

                        // using the estimated constant C, an amplitude can be calculated for the user's tauC
                        // ALISTAIR Modified 18/10/2012:
                        // Now uses interpolated initial values calculated above
                        amplitudeReal = initialRealAmpForUserTauC * Math.Exp(-1 * avgConstantCReal * userTauC * currentTimepoint);
                        amplitudeImag = initialImagAmpForUserTauC * Math.Exp(-1 * avgConstantCImag * userTauC * currentTimepoint);

                        // Generate a new point based on the current time and calculated real and imag amplitudes,
                        // then add the point to the array representing the simulated FID.
                        userTauCInterpolatedPoints[timepointIndex] = new Point(currentTimepoint, amplitudeReal, amplitudeImag);

                    } // end loop through interpolated timepoints

                    // Diagnostic output for debugging
                    /*                    
                    String diagnosticOutput = "User TauC Interpolated Function\n";
                    diagnosticOutput += "TauC," + userTauC + "\n";
                    diagnosticOutput += "Timepoint,Real Value,Imag Value\n";
                    for (int k = 0; k < userTauCInterpolatedPoints.Length; k++)
                    {
                        diagnosticOutput += userTauCInterpolatedPoints[k].T.ToString() +
                                            "," + userTauCInterpolatedPoints[k].V.Real.ToString() +
                                            "," + userTauCInterpolatedPoints[k].V.Imaginary.ToString() + "\n";
                    }
                    System.Diagnostics.Debug.WriteLine(diagnosticOutput);
                    */                  

                    // ALISTAIR: COMMENT ADDED 17/03/2012
                    // The analytics are performed again for the user specified tauC interpolated data.
                    DateTime final_beg = DateTime.Now;
                    this.ExecutionStatus = string.Format("Processing Dimension {0}.", bS);
                    SimulatedData dataForUserTauC = new SimulatedData();
                    dataForUserTauC.TauC = userTauC;

                    // In order to re-use the CalculateBenchmarks() method, the interpolated data needs to be
                    // reconstructed into time/value data in the format imported from the XML data. 
                    for (int timepoint = 0; timepoint < userTauCInterpolatedPoints.Length; timepoint++)
                    {
                        dataForUserTauC.TimeValueRealTable += userTauCInterpolatedPoints[timepoint].T.ToString()
                         + " " + userTauCInterpolatedPoints[timepoint].V.Real.ToString() + "\n";
                        dataForUserTauC.TimeValueImagTable += userTauCInterpolatedPoints[timepoint].T.ToString()
                         + " " + userTauCInterpolatedPoints[timepoint].V.Imaginary.ToString() + "\n";
                    }

                    // add the new tauC data to the list of data for the current indirect dimension
                    bS.DimensionData.Add(dataForUserTauC);

                    // re-run CalculateBenchmarks method for user specified tauC data
                    BenchmarkResult bM_forUserTauC = this.CalculateBenchmarks(bS, dataForUserTauC);

                    // calculate time taken
                    bM_forUserTauC.ExecutionTime = DateTime.Now - final_beg;

                    // insert user specified tauC analysis results to the start of the array of results
                    this.Results.Insert(dimensionNumber-1, bM_forUserTauC);

                } // end foreach BenchmarkSetup
                
            }
            finally
            {
                this.Terminated = true;
            }

        }  // end RunEverything

        // ALISTAIR: COMMENT ADDED 24/01/2011
        // Initiate a new thread for calculation of benchmark statistics.
        // Since these calculations may take a long time, running in a new
        // thread allows the user to continue using the NMRplus website
        // while the calculations are taking place.
        public void Start()
        {
            this.Thread = new Thread(this.RunEverything);
            this.Thread.Priority = ThreadPriority.Lowest;

            this.ExecutionStatus = string.Empty;

            this.Terminated = false;
            this.Thread.Start();
        }

        // ALISTAIR: COMMENT ADDED 24/01/2011
        // Cancels any currently executing benchmark calculations
        // Initiated by clicking "Cancel Benchmarks" button, shown
        // during execution of benchmark calculations
        public void Cancel()
        {
            this.Terminated = true;
            this.Thread.Join();
        }

        // ALISTAIR: COMMENT ADDED 24/01/2011
        // Default constructor
        public BenchmarksBuilder()
        {
            this.ExecutionStatus = string.Empty;
            this.Results = new List<BenchmarkResult>();

            this.Terminated = true;
        }


    }


    // ALISTAIR: 16/08/2011
    // Modified BenchmarkSetup to allow Multi-dimensional data
    // Assign input values to appropriate variables
    [Serializable]
    public class BenchmarkSetup
    {
        /// <summary>
        /// The name of dimension.
        /// </summary>
        [XmlAttribute("Dimension")]
        public string Dimension { get; set; }

        // ALISTAIR Added 04/04/2012
        /// <summary>
        /// A string that indicates how the time domain points
        /// parameter is labeled in the pulse program setup script
        /// </summary>
        //[XmlAttribute("DimPrefix")]
        //public string DimPrefix { get; set; }

        /// <summary>
        /// Number of data points in the time domain for this dimension
        /// </summary>
        [XmlIgnore]
        public int TimeDomain { get; set; }

        /// <summary>
        /// Increment value
        /// </summary>
        [XmlIgnore]
        public double Increment { get; set; }
        
        // ALISTAIR: Added on 26/03/2012
        /// <summary>
        /// Number of data points in the freq domain for this dimension
        /// </summary>
        [XmlIgnore]
        public int NoOfPointsSI { get; set; }

        /// <summary>
        /// User Defined Rotational Correlation Time (TauC)
        /// </summary>
        // Input directly from textbox, not XML tag,
        // so should be ignored by XML reader
        [XmlIgnore]
        public double UserTauC { get; set; }

        /// <summary>
        /// User Defined Signal to Noise from 1D Spectrum
        /// </summary>
        // Input directly from textbox, not XML tag,
        // so should be ignored by XML reader
        [XmlIgnore]
        public double OneDimSigToNoise { get; set; }

        // ALISTAIR: Added on 26/03/2012
        /// <summary>
        /// Number of data points in the freq domain for this dimension
        /// </summary>
        [XmlIgnore]
        public int NoOfTransients { get; set; }

        // ALISTAIR 16/08/2011:
        // Each BenchmarkSetup (Dimension)
        // should contain a list of "Real" datasets
        // and a list of "Imaginary" datasets
        // (one for value of TauC)
        public List<SimulatedData> DimensionData = new List<SimulatedData>();

    }  // end BenchmarkSetup class


    public class SimulatedData
    {

        // ALISTAIR 18/08/2011
        // The value of the Real tag, containing all time/value pairs
        //[XmlIgnore]
        public string TimeValueRealTable { get; set; }

        // ALISTAIR 18/08/2011
        // The value of the Imag tag, containing all time/value pairs
        //[XmlIgnore]
        public string TimeValueImagTable { get; set; }

        // "Real" TauC value
        public double TauC { get; set; }

        // save the maximum time value (which is always the last)
        public double Tmax { get; set; }

        // ALISTAIR: COMMENT ADDED 24/01/2011
        // Extracts all of the transfer function data, within the <BenchmarkSetup> tag of the XML.
        // Each line within the XML file becomes a Point, resulting in an array of Points
        public Point[] GetTransferFunction()
        {
            // ALISTAIR: COMMENT ADDED 24/01/2011
            // Split the input table into an array of strings, delimited by the newline, return or space chars.
            // This will separate x and y values onto consecutive lines.
            // FOR ALL REAL & IMAGINARY SETS
            // NEED TO SPLIT INTO SEPARATE SETS FOR EACH REAL AND IMAGINARY TAG
            string[] tFRealLines = this.TimeValueRealTable.Split(new char[] { '\n', '\r', ' ' }, StringSplitOptions.RemoveEmptyEntries);
            string[] tFImagLines = new string[0];
            // Check if imaginary data exists in the XML before attempting to populate the array
            if (this.TimeValueImagTable == String.Empty)
            {
                tFImagLines = this.TimeValueImagTable.Split(new char[] { '\n', '\r', ' ' }, StringSplitOptions.RemoveEmptyEntries);
            }

            List<Point> points = new List<Point>();

            // Traverse through the array, creating a new Point for each x,y value pair
            // Note that x and y are on alternate lines, so interator is incremented by 2 for each point
            for (int i = 0; i < tFRealLines.Length - 1; i += 2)
            {
                double x, y_real, y_imag;
                x = double.Parse(tFRealLines[i]);  // x becomes the time
                y_real = double.Parse(tFRealLines[i + 1]);  // y_real becomes the real component of complex value
                // check if an imaginary component exists for the corresponding real value
                // if so, assign the current element of the imaginary array to be the imaginary component
                // otherwise, assign a value of zero to the imaginary component
                if (tFImagLines.Length > i+1)
                {
                    y_imag = double.Parse(tFImagLines[i + 1]);  // y_imag becomes the imaginary component of complex value
                }
                else
                {
                    y_imag = 0.0;
                }
                Point p = new Point(x, y_real, y_imag);
                points.Add(p);
            }

            // Assign Tmax as final time in tFLines
            // Final line is final value (Length-1)
            // Second last is final time - Tmax (Length-2)
            Tmax = double.Parse(tFRealLines[tFRealLines.Length - 2]);

            return points.ToArray();

        }  // end GetTransferFunction

    }  // end SimulatedData class


    // ALISTAIR: COMMENT ADDED 24/01/2011
    // Output variables for benchmarking results
    [Serializable]
    public class BenchmarkResult
    {

        public double TauC { get; set; }

        public string Dimension { get; set; }

        public TimeSpan ExecutionTime { get; set; }

        /// <summary>
        /// Number of data points in this dimension
        /// </summary>
        public int TimeDomain { get; set; }

        public PhysicalValue<double> Signal { get; set; }

        public PhysicalValue<double> Noise { get; set; }

        public PhysicalValue<double> Linewidth { get; set; }

        public PhysicalValue<double> SignalToNoise { get; set; }

        public double StNtTD { get; set; }

        public Point[] Spectrum { get; set; }

// ALISTAIR 27/03/2012:
// Temporarily add additional values, for testing purposes
public double StoNfrom1D { get; set; }
public double SI { get; set; }
public double NoOfScans { get; set; }
public double StoNdivRootTD { get; set; }

    }

    [Serializable]
    public struct Point
    {
        public double T { get; set; }

        public Complex32 V { get; set; }

        public Point(double x, double y_real, double y_imag)
            : this()
        {
            this.T = x;
            this.V = Complex32.WithRealImaginary((float)y_real,(float)y_imag);
            //this.V = Complex.FromRealImaginary(y_real, y_imag);
        }
    }

    [Serializable]
    public struct PhysicalValue<T>
    {
        public T Value { get; set; }
        public T Error { get; set; }

        public PhysicalValue(T val, T err)
            : this()
        {
            this.Value = val;
            this.Error = err;
        }
    }

} // end IdeaNMR.Shared namespace
